clear all;

lambda_SOC=0.1*0.5;

spin_Delta=0.0;
spin_nvec=[0,0,1];
spin_nvec=spin_nvec/norm(spin_nvec);

sigx=[0,1;1,0];
sigy=[0,-i;i,0];
sigz=diag([1,-1]);
sig0=eye(2);

gen_SOC_term=@(bdvec) sqrt(2)*(i)*lambda_SOC*(bdvec(1)*sigx+bdvec(2)*sigy+bdvec(3)*sigz);


a1=[2.7508   -1.5882    4.4921];
a2=[0    3.1764    4.4921];
a3=[-2.7508   -1.5882    4.4921];

Rh_pos=[0.125000 0.125000 0.625000 ;
    0.625000 0.125000 0.125000 ;
    0.125000 0.625000 0.125000 ;
    0.125000 0.125000 0.125000];

rec_a1=a1+a2-a3;
rec_a2=a1+a3-a2;
rec_a3=a3+a2-a1;

Rh_pos=Rh_pos-0.125;

Rh_posC=Rh_pos;
Rh_posC(:,1)=Rh_pos(:,1)*a1(1)+Rh_pos(:,2)*a2(1)+Rh_pos(:,3)*a3(1);
Rh_posC(:,2)=Rh_pos(:,1)*a1(2)+Rh_pos(:,2)*a2(2)+Rh_pos(:,3)*a3(2);
Rh_posC(:,3)=Rh_pos(:,1)*a1(3)+Rh_pos(:,2)*a2(3)+Rh_pos(:,3)*a3(3);



Rh_posC


tetra1_posC=Rh_posC
tetra1_C=mean(tetra1_posC,1)

tetra2_posC=-Rh_posC
tetra2_C=mean(tetra2_posC,1)


newx_axis=zeros(4,3);
newy_axis=zeros(4,3);
newz_axis=zeros(4,3);
for ind=1:4
    zvec=tetra1_posC(ind,:)-tetra1_C;
    newz_axis(ind,:)=-zvec/norm(zvec);
end
newz_axis

newx_axis(4,:)=[1,0,0];
newx_axis(3,:)=[-1,0,0];
newx_axis(2,:)=[0.5,0.5*sqrt(3),0];
newx_axis(1,:)=[0.5,-0.5*sqrt(3),0];


for ind=1:4
    xvec=newx_axis(ind,:);
    zvec=newz_axis(ind,:);
    newy_axis(ind,:)=cross(zvec,xvec);
end


dot(newy_axis(2,:),newz_axis(2,:))




tetra1=Rh_posC;
tetra1C=mean(tetra1,1);

tetra2=Rh_posC;
tetra2(1,:)=tetra2(1,:)-a3;
tetra2(2,:)=tetra2(2,:)-a1;
tetra2(3,:)=tetra2(3,:)-a2;

tetra1
tetra1C=mean(tetra1,1)
tetra2
tetra2C=mean(tetra2,1)




figure(1);
clf;


hold on;
linecc= [.5 .5 .5];
plot_tetra_v4(Rh_posC,0*a1,linecc,2)
plot_tetra_v4(Rh_posC,a1,linecc,2)
plot_tetra_v4(Rh_posC,a2,linecc,2)
plot_tetra_v4(Rh_posC,a3,linecc,2)


%plot_tetra(-Rh_posC,0*a1,'r')
plot_tetra_v4(-Rh_posC,a1,linecc,2)
plot_tetra_v4(-Rh_posC,a2,linecc,2)
plot_tetra_v4(-Rh_posC,a3,linecc,2)

plot_tetra_v4(-Rh_posC,a1+a2,linecc,2)
plot_tetra_v4(-Rh_posC,a2+a3,linecc,2)
plot_tetra_v4(-Rh_posC,a1+a3,linecc,2)



hexpts1=[];
hexpts1(1,:)=Rh_posC(3,:)+a1;
hexpts1(2,:)=Rh_posC(2,:)+a2;
hexpts1(3,:)=Rh_posC(1,:)+a2;

hexpts1(4,:)=Rh_posC(3,:)+a3;
hexpts1(5,:)=Rh_posC(2,:)+a3;
hexpts1(6,:)=Rh_posC(1,:)+a1;


hexpts2=[];
hexpts2(1,:)=Rh_posC(4,:)+a1;
hexpts2(2,:)=Rh_posC(1,:)+a1;
hexpts2(3,:)=Rh_posC(2,:)+a3;
hexpts2(4,:)=Rh_posC(4,:)+a3;
hexpts2(5,:)=Rh_posC(1,:);
hexpts2(6,:)=Rh_posC(2,:);



hexpts3=[];
hexpts3(1,:)=Rh_posC(4,:)+a2;
hexpts3(2,:)=Rh_posC(2,:)+a2;
hexpts3(3,:)=Rh_posC(3,:)+a1;
hexpts3(4,:)=Rh_posC(4,:)+a1;
hexpts3(5,:)=Rh_posC(2,:);
hexpts3(6,:)=Rh_posC(3,:);



hexpts4=[];
hexpts4(1,:)=Rh_posC(4,:)+a3;
hexpts4(2,:)=Rh_posC(3,:)+a3;
hexpts4(3,:)=Rh_posC(1,:)+a2;
hexpts4(4,:)=Rh_posC(4,:)+a2;
hexpts4(5,:)=Rh_posC(3,:);
hexpts4(6,:)=Rh_posC(1,:);



sig_x=[0,1;1,0];
sig_y=[0,-i;i,0];
sig_z=[1,0;0,-1];
sig_0=eye(2);

gen_spinh=@(nnvec) (sig_x*nnvec(1)+sig_y*nnvec(2)+sig_z*nnvec(3))/norm(nnvec);

eva_spinS =@(state_wf) [state_wf'*sig_x*state_wf,state_wf'*sig_y*state_wf,state_wf'*sig_z*state_wf];



newnvec_tetra=Rh_posC(3,:)-tetra1C;
newnvec_tetra=newnvec_tetra/norm(newnvec_tetra);





tripts=Rh_posC(1:3,1:3);
tripts_xaxis=newx_axis(1:3,1:3);

frame_x=newx_axis(1:3,1:3);
frame_y=newy_axis(1:3,1:3);


%quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    plt_newD2r_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% 
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end

for ind=1:3
    [dot(tripts_xaxis(ind,:),frame_x(ind,:)),dot(tripts_xaxis(ind,:),frame_y(ind,:))]
end





shiftvec=a2;
tripts=Rh_posC([1,2,4],1:3);
%tripts_xaxis=newx_axis([1,2,4],1:3);

newnvec=Rh_posC(3,:)-tetra1_C;
newnvec=newnvec/norm(newnvec);
newnvecx=[1,0,0];
newnvecy=cross(newnvec,newnvecx);
tripts_xaxis=zeros(3,3);
tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
tripts_xaxis(3,:)=[-1,0,0];

tripts_xaxis=-tripts_xaxis;

frame_x=newx_axis([1,2,4],1:3);
frame_y=newy_axis([1,2,4],1:3);

for ind=1:3
    plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
%     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
%     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end

for ind=1:3
    %[dot(tripts_xaxis(ind,:),frame_x(ind,:)),dot(tripts_xaxis(ind,:),frame_y(ind,:))]
end


%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')


R3mat=zeros(3);
R3mat(3,3)=1;
R3mat(1,1)=-0.5;
R3mat(2,2)=-0.5;
R3mat(1,2)=0.5*sqrt(3);
R3mat(2,1)=-0.5*sqrt(3);


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a1;
tripts=Rh_posC([3,1,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')

frame_x=newx_axis([3,1,4],1:3);
frame_y=newy_axis([3,1,4],1:3);

for ind=1:3
    plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
end


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a3;
tripts=Rh_posC([2,3,4],1:3);

frame_x=newx_axis([2,3,4],1:3);
frame_y=newy_axis([2,3,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
end




linecolor=[0.0,0.4,0.0];
linewidth=7;

tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:);
tri_pts(2,:)=Rh_posC(2,:);
tri_pts(3,:)=Rh_posC(3,:);
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;

plot_polygon_line(tri_pts,linecolor,linewidth)

tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:)+a1;
tri_pts(2,:)=Rh_posC(3,:)+a1;
tri_pts(3,:)=Rh_posC(4,:)+a1;
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
plot_polygon_line(tri_pts,linecolor,linewidth)


tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:)+a2;
tri_pts(2,:)=Rh_posC(2,:)+a2;
tri_pts(3,:)=Rh_posC(4,:)+a2;
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
plot_polygon_line(tri_pts,linecolor,linewidth)


tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(2,:)+a3;
tri_pts(2,:)=Rh_posC(3,:)+a3;
tri_pts(3,:)=Rh_posC(4,:)+a3;
plot_polygon_line(tri_pts,linecolor,linewidth)


linecolor=[0.0,0.8,0.0];
linecolor='cyan';
linecolor='#FF8800';
 linecolor=[0.9290, 0.6940, 0.1250];
 linecolor=[0.75, 0, 0.75]; % purple
  linecolor=='m';
  linecolor='cyan';

pt1=Rh_posC(2,:);
pt2=Rh_posC(4,:)+a1;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
pt1=Rh_posC(3,:);
pt2=Rh_posC(4,:)+a2;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(1,:);
pt2=Rh_posC(4,:)+a3;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(3,:)+a1;
pt2=Rh_posC(2,:)+a2;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(1,:)+a2;
pt2=Rh_posC(3,:)+a3;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)


pt1=Rh_posC(2,:)+a3;
pt2=Rh_posC(1,:)+a1;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)




coorcenter=[-4.5,-3.,1];
coordinatelen=1.5;
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% %quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)


axisvec = rec_a1/norm(rec_a1)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)
axisvec = rec_a2/norm(rec_a2)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)
axisvec = rec_a3/norm(rec_a3)*coordinatelen;
%quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',15,'Color','k','LineWidth',7)





hold off;


box on;
axis equal;


%view([0,0,1])
view([0,13])

set(gca,'visible','off')
material dull

 camlight('headlight')






%%


set(gca,'visible','off')
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4h-D2CLS.eps','ContentType','image','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'





%%


figure(2);

clf;

hold on;


coorcenter=[-4.5,-3.,1];
coordinatelen=1.5;
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% %quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)

axisvec = rec_a1/norm(rec_a1)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)
axisvec = rec_a2/norm(rec_a2)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)
axisvec = rec_a3/norm(rec_a3)*coordinatelen;
quiver3(coorcenter(1),coorcenter(2),coorcenter(3),axisvec(1),axisvec(2),axisvec(3),'MaxHeadSize',40,'Color','k','LineWidth',15)



set(gca,'visible','off')

view([0,13])


material dull
 camlight('headlight')

hold off;


box on;
axis equal;

%%




set(gca,'visible','off')
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4h-D2CLS-hex-coordinate-frame.eps','ContentType','vector','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'






